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ABSTRACT 

In the past five years, approximately one third of the 65 pulsars discovered by radio observations of 
Fermi unassociated sources are black widow pulsars (BWPs). BWPs are binary millisecond pulsars 
with companion masses ranging from 0.01-0.1 solar masses which often exhibit radio eclipses. The 
bloated companions in BWP systems exert small torques on the system causing the orbit to change on 
small but measurable time scales. Because adding parameters to a timing model reduces sensitivity to 
a gravitational wave (GW) signal, the need to fit many orbital frequency derivatives to the timing data 
is potentially problematic for using BWPs to detect GWs with pulsar timing arrays. Using simulated 
data with up to four orbital frequency derivatives, we show that fitting for orbital frequency derivatives 
absorbs less than 5% of the low frequency spectrum expected from a stochastic gravitational wave 
background signal. Furthermore, this result does not change with orbital period. Therefore, we 
suggest that if timing systematics can be accounted for by modeling orbital frequency derivatives and 
is not caused by spin frequency noise, pulsar timing array experiments should include BWPs in their 
arrays. 


1. INTRODUCTION 

Black widow pulsars (BWPs) are a special class of 
binary millisecond pulsars (MSPs). The y often have 
orbit al periods ranging from 2-20 hours (Chen et al. 
2013). However, the defining feature of BWPs is their 
low mass companions, with masses much less than 0.1 


Mq ( [Roberts pQ13 ). 

Even though BwPs have very small companions, the 
i onized gas from the co mpanions often eclipses the pulsar 


(Roberts et al. 
wind blows stel 


2014). This occurs because the pulsar 
ar material away from the companion. 


but leaves plasma in th e system, creating a screen for the 
pulsar’s radio emission ( Thompson [1995 ). The eclipse ir 
the first discovered black widow s ystem, rSR B1957-1-20 


covers about 10% o f the orbit (Fruchter et al. 
Roberts et al. 12^ . 


1988 


However, in this work, we are most concerned with 
the stochastic variation in orbital period exhibited by 
BWP systems. In fact, in a few years of observing 
PSR B1957+ 20, the sign of the orbital period deriva¬ 


tive changed (Arzoumanian et al. 1994). In order to 


explain this observation, Applegate bnaham ([1994 ) 
gave a mechanism for the changes in orbital period over 
time. Their model starts from the fact that there are 
large tidal forces between the pulsar and companion due 
to their tight orbit. These tidal forces drive stellar con¬ 
vection in the companion. Gombined with rapid rota¬ 
tion, stellar convection gives rise to a strong magnetic 
field. This strong magnetic field and differential rota¬ 
tion in the compan ion cause the com panion to become 
more or less oblate ( Applegat^|1992 ). If the companion 
becomes more oblate, the centrepital force between the 
pulsar and the companion increases while their angular 
momentum remains constant, causing a decrease in or¬ 
bital period. If the companion becomes less oblate, the 
centrepital force decreases causing the orbital period to 


increase. 

Due to the high orbital variability of BWP systems, 
precisely modeling the pulse arrival times for the pulsar 
often requires modeling the orbit with many orbital fre¬ 
quency derivatives (OFDs) in our timing solutions. How¬ 
ever, fitting many OFDs reduces sensitivity to gravita¬ 
tional waves (GWs). Because the parameters of a pulsar 
are determined from the same data used to search for 
GWs, some of the GW signal will be absorbed in fits for 
pulsar timing parameters, as the GW signal will appear 
in at leas t a small way to be cov ariant with the fit pa¬ 


rameters (Demorest et al. 2013). This process reduces 


the amount of GW signal in the data after the timing fit 
r educing sensitivity t o a GW source. 


Nice et al. (2000) showed PSR B1957-F20, the first 


discovered BWP system, had what appeared to be signif¬ 
icant stochastic deviation from the best fit timing model, 
or timing noise, over many years and suggested that the 
source of this noise is intrinsic to the system, such as 
stochastic variation in spin period. This result has led 
pulsar timing array (PTA) collaborations to be reluctant 
about including BWPs in PTAs to search for GWs. The 
variations in orbital period make BWPs difficult to time 
and perhaps contribute to the timing noise. However, if 
we can account for the timing variations in a way that 
does not significantly reduce sensitivity to GWs, then 
black widow pulsars should be considered for inclusion in 
PTAs. We show that modeling orbital frequency deriva¬ 
tives does not significantly reduce sensitivity to GWs and 
suggest that other timing variations can be explained by 
variations in dispersion measure. 

All of this comes at a time when BWPs are being 
discovered in droves. One of the most successful ways 
of searching for MSPs has be en surveying y-ray sources 


detect ed by the Fermi-LAT (Ray et al 


2012 Roberts 


2013). Among these newly discovered 4ernii MSPs 





































2 


about a third of th e 65 pulsars discovered usin^ thi s 


method are BWPs (Ray et ah 2012 [Roberts 2013) 


Given that BWPs are being discovered at a very high 
rate, it makes sense to revisit the question of whether 
BWPs can be timed precisely enough to use in PTAs. 

2. METHODS 

2.1. Quantifying the Amount of Gravitational Wave 
Signal Post-fit 

When a GW passes between the Earth and a pulsar, 
the distance between the pulsar and th e Earth is altered 
as a function of time (Detweiler ||1979). Therefore, there 
will b e timing variations from GWs. Jaffe & Backer 
(2003) showed that the characteristic strain spectrum for 


a stochastic background of supermassive black hole bina¬ 
ries will be a power law of index -2/3. Thus, we can esti¬ 
mate the effect of such a background on a pulsar’s times 
of arrival (TOA). By hypothesizing that all of the struc¬ 
ture in the timing residuals is due to gravitational waves, 
we c an place a liniit on the amplitude of this background 
(e.g. Kaspi et al. 

In order to quantify now much gravitational wave sig¬ 
nal is absorbed by fittin g orbital frequency deriva tives, 


we use a method from Demorest et al. (2013) that 


assumes there is a stochastic gravitational wave back¬ 
ground signal of spectral index -2/3 in a data set and cal¬ 
culates how much of that assumed signal will be present 
in the post-fit timing residuals. Eirst, we calculate the 
pre-fit gravitational wave covariance matrix , (see 

[Demorest et al. ||2Q13 van Haasteren et al. 2QQ9[ ). This 


is done by calculating the statistical expectation value 
of ya{U)ya{tj) for all pulse TO As i and j, where ya{t) 
is the amount a gravitational wave background signal 
shifts a pulse TOA. Then, we apply the timing fit to 
this matrix to get the post-fit, or residual, gravitational 
wave covariance matrix, Next, we weight 

by the diagonal matrix of pulse TOA uncertainties, W, 
to get WCp^W . Then, we diagonalize WCp^W to 
form an orthonormal basis of eigenvectors that is com¬ 
pletely orthogonal to the timing fit, optimally capturing 
the gravitational wave signal left in the residuals after 
the fit. 

The eigenvalues Xi are the square of the amplitude each 
element of this basis contributes to the timing residuals. 
By summing over y/Xi^ we compute the amplitude a GW 
background signal in our data would contribute to the 
timing residuals. By repeating this process while fitting 
only for spin frequency and spin frequency derivative, we 
calculate the total amount of GW background signal ef¬ 
fectively in our data. Any timing fit will include spin 
frequency and spin frequency derivative, so the signal 
lost due to fitting those two parameters is effectively not 
in our data. This is shown in Equation 1, where n is the 
number of TOAs, Af®® are the eigenvalues of WCp^W 
where all timing parameters were free, and are the 

eigenvalues of WCp^W where all the timing parame¬ 
ters were fixed except the spin frequency and spin fre¬ 
quency derivative. %GWs represents the percentage of 
assumed GWs present in the data after the fit. 


%GWs = 




2.2. Sensitivity vs. Orbital Period 

In order to investigate how sensitivity to GWs c hanges 
with orbital period, using the TEMP02 plugin FAKE (Hobbs 


et al. |2006 ), we simulated ten data sets with timing vari- 
ations of rms 100 ns with no orbital frequency deriva¬ 
tives from PSR J0023-F0923. It should be noted that 
PSR J0023+0293 shows very little orbital period v aria- 


zoumanian et al. [| 

2015 

Roberts et al. [ 

12014). Eachc 

.at a 

set ranged from 5^ 

1000 MJD to 54000 iV 

IJD. We repeated 


X 100% 


( 1 ) 


this process, changing the orbital period within the range 
2.7 hours to 300 years. We then fit each set of simulated 
TOAs to its respective model with the correct orbital 
period. Eor each fit, we calculated the percent of a GW 
background signal left in the data after the fit and took 
the average value of the ten data sets with the same or¬ 
bital period as the correct value. We then repeated this 
process fitting up to ten orbital frequency derivatives. 
The results are shown in Eigure 1. 

2.3. Sensitivity vs. Number of Modeled Orbital 
Frequeney Derivatives 

We used simulated data to quantify the effect on GW 
sensitivity of modeling many orbital frequency deriva¬ 
tives. We simulated five data sets from PSR J1748- 
2446ad, PSR J1748-2021D, PSR J1748-2446P, and PSR 
J2129-04 with rms 10“^ from MJD 51000 to 54000. We 
chose these pulsars because each of them has a timing so¬ 
lution and requires modeling of several orbital frequency 
derivative terms. Eor each pulsar, each data set has a 
different number of orbital frequency derivative signals, 
ranging from no derivatives to four derivatives. We then 
made 17 model files for each data set, each modeling 
the simulated data with a different number of orbital 
frequency derivatives ranging from 0 to 16 derivatives. 
This may sometimes over-fit our simulated observations. 
We fit each model to each data set and computed the 
amount of GW background signal in the data after the 
fit for each fit. The results are shown in Eigure 2. 

2.4. Transmission Speetrum of PSR J1748-244dP 

In order to further demonstrate that the signal re¬ 
moved from modeling orbital frequency derivatives looks 
nothing like the signal of the GWs PTAs might see, we 
calculate the transmission spectrum of GWs for a tim¬ 
ing model of based off PSR J1748-2446P, but with up to 
four OEDs modeled. Eirst, we simulated pulse TOAs at 
different dates ranging from 51000 MJD to 54000 MJD, 
such that on average there is an observation every two 
weeks with white timing noise of rms 100 ns with up 
to four OED signals in the data. Then, for each of 
the four timing models of PSR J1748-2446P with be¬ 
tween one and four OEDs modeled, we fit the simulated 
TOAs to each corresponding timing model with the cor¬ 
rect number of OEDs modeled. We then made a nor¬ 
malized Lomb-Scargle periodogram from the residuals of 
each timing fit. This process was repeated for 100 differ¬ 
ent sets of simulated TOAs in order to minimize noise. 
We plot the average Lomb-Scargle periodogram for each 
number of OEDs fit in Eigure 3. This shou ld correspond 
to the tra nsmission spectrum described in[Blandford et 

ar[(fT^. 


3. RESULTS 
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Spin 

Orbital 


%GWs in 

Residuals 







Pulsar 

Period 

(ms) 

Period 

(hrs) 

Observed 

OFDs 

0 OFD signals 

0 OFDs fit 

4 OFD signals 
16 OFDs fit 

% Ghange 

PSR J1748-2021D 

13.50 

6.87 

2 

96.6 

93.9 

2.7 

PSR J1748-2446ad 

1.40 

26.23 

6 

97.9 

94.6 

3.4 

PSR J1748-2446P 

1.73 

8.70 

11 

97.9 

94.9 

3.1 

PSR J2129-04 

7.61 

15.25 

3 

99.8 

94.1 

5.7 

Average 






3.7 


TABLE 1 

The percent change in the amount of GWs that survive the fit after adding four orbital frequency derivative signals 

AND FITTING OUT 16 ORBITAL FREQUENCY DERIVATIVE TERMS. 


3.1. Sensitivity vs. Orbital Period 

We found that there are three GW sensitivity regimes 
in orbital period: a high sensitivity regime, a low sensi¬ 
tivity regime, and a transition regime. The low sensitiv¬ 
ity regime covers orbital periods greater than the length 
of the data set, the transition regime covers orbital peri¬ 
ods slightly less than the length of the data set, and the 
high sensitivity regime includes all orbital periods signif¬ 
icantly shorter than the data span. This is expected as 
when the orbital period is greater than the length of the 
timing data set, GWs and orbital parameters are com¬ 
pletely covariant. 

As shown in Figure 1, in the low sensitivity regime 
at the length of the data set, there is rougly a 65-75% 
decrease in sensitivity between modeling no orbital fre¬ 
quency derivatives and modeling 10 orbital frequency 
derivatives. In the transition regime, there is roughly a 
50-60% decrease in sensitivity between modeling no or¬ 
bital frequency derivatives and modeling 10 orbital fre¬ 
quency derivatives. In the high sensitivity regime, there 
is at most a 10% decrease in sensitivity between mod¬ 
eling no orbital frequency derivatives and modeling 10 
orbital frequency derivatives. 
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Fig. 2.— Amount of a GW background signal not absorbed post- 
fit vs. number of orbital frequency derivatives fit for data sets with 
different numbers of orbital frequency derivative signals. We used 
simulated data from a timing solution of PSR J1748-2446ad (or¬ 
bital period = 26.23 hrs) with the respective number of orbital 
frequency derivative signals. 


3.2. Sensitivity vs. Number of Modeled Orbital 
Frequeney Derivatives 

Figure 2 shows that modeling orbital frequency deriva¬ 
tives has very little effect on sensitivity to GWs and 



Fig. 1. — Amount of a GW background signal remaining post-fit 
vs. orbital period for different numbers of orbital frequency deriva¬ 
tive terms fit. We used simulated data from a timing solution for 
PSR J0023+0293 without orbital frequency derivative terms. 


that the decrease in sensitivity due to modeling orbital 
frequency derivatives does not strongly depend on how 
many orbital frequency derivative signals are in the data. 
For PSR J1748-2446ad, there is a decrease in GW sensi¬ 
tivity of only 3.4% between fitting no orbital frequency 
derivatives to a data set with no orbital frequency deriva¬ 
tives and fitting 16 orbital frequency derivatives to a data 
set with four orbital frequency derivatives. Furthermore, 
when this analysis was repeated for other pulsars, the 
decrease in GW sensitivity varied only slightly, giving an 
average decrease of 3.7%. None of the pulsars showed a 
significant decrease in GW sensitivity due to modeling 
more orbital frequency derivatives. 


4. DISGUSSION 

4.1. Sensitivity vs. Orbital Frequeney 

We expect to lose sensitivity to GWs when the orbital 
period is greater than the length of the data set. Because 
PT As are most sens itive to frequencies on the order of 
^ ( Ellis et aL~||2Q12 ), where T is the length of the data 
set, the full period of any GW detected will not be in 
the data. Therefore, when we search our data for GWs, 
we are searching in part for a signal that is an incom¬ 
plete sine wave. As the length of the orbit approaches 
the length of the data set, the situation is analogous to 
the dramatic decrease in GW sensitivity when fitting 
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Fig. 3.— GW Transmission spectrum for PSR J1748-2446P. Both panels show the normalized Lomb-Scargle periodogram of the timing 
fit residuals for simulated data from PSR J1748-2446P with up to four OFD signals. The left panel shows the entire spectrum while the 
right panel is zoomed in on the orbital frequency. Each color corresponds to the number of OFDs fit for in the model. The dashed vertical 
lines correspond to the inverse of the length of the data set, the frequency corresponding to a period of a year, and the orbital frequency. 
It is clear that as more OFDs are fit, the dip in Fourier power around the orbital frequency is spread out around the orbital frequency. 


for many spin period derivatives. The timing data set 
will not contain a full orbit, meaning the signal due to 
the orbit will be indistinguishable from the signal due to 
GWs. Thus, the GW signal and the signal from the or¬ 
bit become covariant. In other words, any GW signal in 
our data will appear to be an error in orbital frequency, 
meaning the fit will absorb the timing variations due to 
GWs into the orbital frequency parameter, greatly re¬ 
ducing sensitivity to GWs. However, orbital periods of 
BWPs range fro m 2-20 hours and o ur timing data sets 
cover many years (Ghen et al. Therefore, we only 

need to be concerned about the high sensitivity regime of 
orbital period, which shows modeling many orbital fre¬ 
quency derivatives has little effect on sensitivity to GWs. 


4.2. Sensitivity vs. Number of Modeled Orbital 
Frequeney Derivatives 

We should not be surprised that modeling orbital fre¬ 
quency derivatives does not significantly reduce sensitiv¬ 
ity to GWs. The signal due to a binary orbit is a sine 
wave. The Fourier transform of a sine wave is a sharp 
spike at the frequency of the sine wave. Even when the 
frequency is slowly changing, as in the orbital frequency 
derivatives, all of the Fourier power remains around the 
orbital frequency. It takes less than a day for the a BWP 
system to complete an orbit, while it takes many years for 
a full GW to pass through the solar system, translating 
to a difference of more than three orders of magnitude 
in period. In other words, because the orbital timescales 
of BWPs and the timescales of GWs are dramatically 
different, the signal due to a binary orbit with changing 
frequency is very different from the signal due to GWs 
that PTAs are attempting to detect. 

We can see this effect in Figure 3, which corresponds to 
the transmission spectrum of GWs for simulated TO As 
from a modified model of PSR J1748-2446P. In this fig¬ 
ure, we see both the transmission spectrum over a wide 
range of frequencies (left panel) and over a narrow win¬ 
dow around the orbital frequency (right panel). The left 
panel shows three large decreases in GW transmission: 
one at frequencies lower than one at 1 year“^, and 
one at the orbital frequency. The decrease in GW trans¬ 


mission at frequencies lower than ^ is expected since 
PTAs are not sensitive to GWs with periods longer than 
the length of the data set. Any GW of with a period 
of one year will be modeled as an error in the pulsar’s 
position, leading to the transmission dip at 1 year“^. 
Similarly, any GW of period equal to the orbital period 
will be modeled as an error in orbital period, producing 
the third decrease in GW transmission at the orbital fre¬ 
quency. The left panel also shows that the orbital period 
is orders of magnitude in frequency away from frequen¬ 
cies of detectable GWs. Furthermore, the right panel 
shows that as more OFDs are modeled, all the Fourier 
power removed is still around the orbital frequency and 
is merely broadened around the orbital frequency. Thus, 
the signals due to OFDs should not be covariant with 
detectable GWs. Because of this, we should not expect 
a large loss in GW sensitivity when accounting for more 
orbital frequency derivatives. Therefore, if the timing 
noise present in BWPs can be accounted for by mod¬ 
eling orbital frequency derivatives, they can be used in 
PTAs. 


4.3. Previous Results and Limitations 


Even though modeling timing noise as orbital fre¬ 
quency derivatives does not remove much sensitivity to 
gravitational waves, the long term timing result s of the 


original B WP showed significant timing noise (Nice et 
al. (|2000|). The largest amplitudes of the timing re sidu¬ 
als were 3' 


0 Rs at 430 MHz and 20 /r s at 575 MHz ([Nice 
et al. 2000). For this timing noise, Nice et al. (2000p 
otters the explanation of dispersion measure (DM) vari¬ 
ations, but suggest that the timing noise lies within the 
system itself. We submit that this timing noise is rea¬ 
sonably explained solely by DM variations. Using the 
equation t = timing delay due to dis¬ 

persion, w here t is the timing delay a nd / the observing 
frequency (Lorimer & Kramer ||2004), we calculate that 
in order to explain the residuals of 30 /is at 430 MHz 
and 20 /rs at 575 MHz, ADMs of -0.0013 pc/cm^ and 


0.0016 pc/cm^ are required, respectively. Furthermore, 
accounting for the timing noise between these two resid¬ 
uals requires the most extreme rate of change in DM. 
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The necessary ADM is 0.0029 pc/cm^ over a period of 
two and a half years, requiring a rate of change in DM of 
0.0012 pc/cm^/yr. This value is well withi n the range of 


published DM variations for P TA MSPs (jArzoumanian 
et al. [K^ et al. |2013] ) 


5. CONCLUSIONS 

Given that the observed long term timing variations 
in PSR B1957+20 can reasonably be explained by DM 
variations, we recommend measuring the DM variations 
in this BWP in order to confirm our hypothesis. Further¬ 
more, in order to test whether or not the timing varia¬ 


tions seen by Nice et al. (2000) were due to DM varia¬ 
tions, we recommend long term multi-frequency timing 
observations of suitable BWPs. If the observed timing 
noise in black widow pulsars can be accounted for with 
DM variations and orbital frequency derivatives, we rec¬ 
ommend using suitable BWPs (i.e. non-eclipsing strong, 
fast, and narrow pulsars) in PTAs. 

BWPs are difficult to time accurately due to their or¬ 
bital dynamics. However, this can be modeled effectively 
by orbital frequency derivatives. Even though we expect 
some GW sensitivity to be lost due to fitting for orbital 
frequency derivatives, we have shown that only a small 
amount of GW sensitivity is lost. 
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